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Dynamic simulation of a beta-type Stirling engine with cam-drive mechanism used in concentrating 
solar power system has been performed. A dynamic model of the mechanism is developed and then 
incorporated with the thermodynamic model so as to predict the transient behavior of the engine in the 
hot-start period. In this study, the engine is started from an initial rotational speed. The torques exerted 
by the flywheel of the engine at any time instant can be calculated by the dynamic model as long as the 
gas pressures in the chambers, the mass inertia, the friction force, and the external load have been 
evaluated. The instantaneous rotation speed of the engine is then determined by integration of the 
equation of rotational motion with respect to time, which in return affects the instantaneous variations 
in pressure and other thermodynamic properties of the gas inside the chambers. Therefore, the transient 
variations in gas properties inside the engine chambers and the dynamic behavior of the engine 
mechanism should be handled simultaneously via the coupling of the thermodynamic and dynamic 
models. An extensive parametric study of the effects of different operating and geometrical parameters 
has been performed, and results regarding the effects of mass moment of inertia of the flywheel, initial 
rotational speed, initial charged pressure, heat source temperature, phase angle, gap size, displacer 
length, and piston stroke on the engine transient behavior are investigated. 

© 2010 Elsevier Ltd. All rights reserved. 


1. Introduction 

Among the most related solar energy technologies, concen¬ 
trating solar power systems with Stirling engines (CSP-Stirling) 
have been treated as one of the most promising energy sources in 
recent years. The CSP-Stirling systems require solar absorbers to 
receive solar radiation and then use the Stirling engines to convert 
the received solar energy to electrical power. Theoretically, the 
performance of the systems for the generation of electrical power is 
dependent on the receiver temperature and the thermal efficiency 
of the engines. It has been demonstrated that the highest system 
efficiency of the CSP-Stirling system reaches 30% in converting 
direct-normal incident solar radiation into electricity after 
accounting the power losses [1], The SES 25-kW SunCatcher [2] 
consisting of a solar dish concentrator, a solar tracker, and a Stir¬ 
ling engine has received great attention from the energy-related 
researchers in recently. In the work of Reinalter et al. [3], the 
performance of a 10-1<W dish/Stirling system was evaluated, which 
was considered to be the latest development step of the EuroDish 
system by improving many components. An axisymmetric 
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computational fluid dynamics approach to the analysis of the 
working process of a solar Stirling engine was performed by 
Mahkamov [4]. Liao [5] also conducted a simulation of the CSP 
system and assessed its feasibility in the Taiwan area. Most recently, 
a number of new articles [6—8] have come to attention regarding 
the CSP technology. 

In the past several decades, numerous researches [9—13] have 
been performed for the development of the Stirling engines. 
A detailed literature review for the solar-powered Stirling engines 
was given by Kongtragool and Wongwises [14]. 

Modeling of the physical transport phenomena in the Stirling 
engines has also received great attention in the past several 
decades. The Schmidt theory [15] has been used widely due to its 
mathematical simplicity. The non-isothermal analysis was first 
carried out by Finkelstein in 1960s [16]. The author evaluated the 
heat transfer coefficient and obtained the heat transfer and the 
temperature variation in the working space. This work was notified 
as one of the major improvement in the modeling of Stirling engine 
processes. Lately, more analytical and simulation models have been 
proposed [17—22] for predicting the performance of the Stirling 
engines. 

By using a quasi-steady flow model based on the works of Urieli 
and Berchowitz [20], thermodynamic analysis of a gamma-type 
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Nomenclature 

V 

Volume (m 3 ) 



W 

Weight of the components (Kg) 

a 

Acceleration (m/s 2 ) 

Wj 

Weight of link 2—3 (Kg) 

Cv 

Constant volume specific heat (kj/kg I<) 

W 2 

Weight of link 5—6 (Kg) 

Cp 

Constant pressure specific heat (kj/kg K) 

W 

Cycle-average output power (W) 

dt 

Time step (s) 

Y 

Displacement (m) 

G 

Regenerative channel gap (m) 



h 

Enthalpy (kj/kg) 

Greek symbols 

I 

Mass moment of inertia (kg m/s 2 ) 

h P 2 

Angle between the central axis and the linkages 2—3 

la 

Length of displacer (m) 


and 5—6 (rad) 

Ip 

Length of piston (m) 

V 

Thermal efficiency 

k 

Height of the engine (m) 

P 

Dynamic viscosity (kg/m s) 

h, h 

Lengths of the link 2—3 and 5—6 of the cam-drive 

6 

Crank angle (rad) 


mechanism (m) 

P 

Fluid density (kg/m 3 ) 

/dt 

Length from the joint 5 to the lower surface of 

(0 

Rotational speed (rad/s) 


displacer (m) 

a 

Angular acceleration (rad/s 2 ) 

m 

Mass (kg) 

£ 

Rod’s angular acceleration (rad/s 2 ) 

m 

Mass flow rate (kg/s) 

£eff 

Regeneration effectiveness 

p 

Pressure (kPa) 

$ 

Phase angle (degree) 

Q 

Volumetric flow rate (m 3 /s) 




Radius of displacer (m) 

Subscripts 

r 2 

Core radius of cylinder (m) 

ave 

Cycle-average 

R 

Gas constant (J/kg K) 

c 

Compression chamber 

Ra 

Offset distance from the displacer crank to the center 

d 

Displacer 


of gear (m) 

e 

Expansion chamber 

R P 

Offset distance from the piston crank to the center of 

eff 

Effectiveness 


gear (m) 

fr 

Friction 

Ra 

Thermal resistance of heating head (°C/1<W) 

i 

Initial value 

Rt 2 

Thermal resistance of cooling jacket (°C/1<W) 

load 

External load 

t 

Time (s) 

P 

Piston 

f p 

Period (s) 

r 

Regenerative channel 

T 

Temperature (K) 

s 

Steady state 

T h 

Heat source temperature (K) 

tin 

Tangential force 

t l 

Heat sink temperature (1<) 

2-3 

link 2—3 

U 

Internal energy (kj/s) 

5-6 

link 5—6 

V 

Velocity (m/s) 




Stirling engine in non-ideal adiabatic conditions was presented by 
Parlak et al. [21], In the work of Karabulut et al. [22], a beta-type 
Stirling engine was built and tested. The beta-type Stirling engine 
was considered to have the largest specific cyclic work as compared 
to other types of the Stirling engines. Most recently, Cheng and 
Yu [23] presented a numerical model for a steadily-operated beta- 
type Stirling engine with rhombic-drive mechanism. By taking into 
account the non-isothermal effects, the effectiveness of the 
regenerative channel, and the thermal resistance of the heating 
head, the authors carried out the parametric study of the depen¬ 
dence of the output power and thermal efficiency on the geomet¬ 
rical and physical parameters, involving regenerative gap, distance 
between two gears, offset distance from the crank to the center of 
gear, and the heat source temperature. 

However, in practices, the study of the transient operating 
process of the engines is essential and necessary, and hence 
dynamic simulation of the engine forms another critical issue of the 
engine’s operation performance. Unfortunately, to the authors’ 
knowledge, most existing simulation model for the beta-type 
Stirling engines analysis only dealt with the steady-state perfor¬ 
mance of the engines, and the transient behavior of the engine in 
the start period has not been sufficiently investigated. 

For an internal combustion engine, Giakoumis and Rakopoulos 
[24] provided a numerical model accounting for the transient 
conditions with the slider-crank-driving mechanism. Lately, the 
same group of authors [25] extended the study to include the 


instantaneous torsional deformation in both the steady-state and 
the transient conditions. The crankshaft angular momentum 
equilibrium was formulated by taking into account instantaneous 
gas, inertia and friction loads as well as the stiffness and damping 
torque contributions. However, the model for the internal 
combustion engines developed in Refs. [24,25] cannot be directly 
applied for the beta-type Stirling engines. 

Under these circumstances, to predict the transient behavior of 
the beta-type Stirling engines with cam-drive mechanism during 
the hot-start period, in the present study a dynamic model is 
developed and then incorporated with the thermodynamic model 
developed by Cheng and Yu in [23] so as to investigate the transient 
variations in gas properties inside the engine chambers and the 
dynamic behavior of the engine mechanism. By combining the 
thermodynamic model [23] and the present dynamic model, an 
extensive parametric study of the effects of mass moment of inertia 
of the flywheel, initial rotational speed, initial charged pressure, 
heat source temperature, phase angle, gap size, displacer length, 
and piston stroke on the engine’s transient behavior is performed, 
and results coming out will be discussed in the subsequent sections. 

2. Dynamic model 

The schematic diagram of the beta-type Stirling engines with 
cam-drive mechanism is illustrated in Fig. 1. The concentric dis¬ 
placer and the piston are placed in one single cylinder. A connecting 
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Fig. 1. Schematic diagram of the beta-type Stirling engine with cam-drive mechanism. 

rod connects the displacer and the flywheel. The joint between the 
rod and the flywheel is distant from the flywheel center at an offset 
distance, Ra. Another rod connects the piston and the flywheel, and 
the joint between the rod and the flywheel is of an offset distance, 
R p . from the flywheel center. The wall of the expansion chamber is 
maintained at a higher temperature (Th) and that of the compres¬ 
sion chamber is at a lower temperature (T l ). The displacer of the 
engine is designed to move the gas medium back-and-forth 
between the two chambers. With the cam-drive mechanism, the 
displacer is able to move with a phase angle ($) ahead of the piston 
by adjusting the positions of the joints connecting the rods and the 
flywheel. The piston delivers the output power of the engine by 
flywheel rotation. All connecting rods of the engine are assumed to 
be rigid bodies with the centers of gravity located in the centers of 
the rods. All the engine components may experiences both rota¬ 
tional and translational movement while the engine is operating. 

The displacement of the piston and displacer with the cam-drive 
mechanism are described as: 

Y p = lip + Rpcosg - 6 ) + licosd?!) (1) 


Yd = Id + fdt + KdCOs(0-'| + $)+J 2 cos(02) (2) 

respectively, where 

0i = sin -1 (^sin (|-0)^ (3) 

and 

02 = sin' 1 (I'sin («-£+*)) (4) 

The reciprocating velocities of the piston and displacer are 
obtained by differentiating Eqs. (1) and (2) with respect to time as 

v p = R p wsin(-— 0) - l 1 w 2 -3sin(0 1 ) (5) 

and 

"d = -fid wsin (0 ~2 + ®)~ ; 2 w 5-6sin(0 2 ) (6) 


respectively, where w (=d0/dt) is the angular speed of the flywheel; 
w 2 -3 (=d|8]/dt) is the angular speed of link 2—3; and 104.5 (=d/? 2 /dt) 
is the angular speed of link 5—6. 


and 
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Furthermore, the accelerations of piston and displacer could be 
obtained by differentiating the velocity equations, and one has 

dp =R p asin^ — ffj - R p m 2 cos- O'j -l^e 2 _ 3 sin(| 8 1 ) 

-/iwl- 3 cos(/?i) (7) 

a d = -^d asin (^-^ + (,) ) -Rd w2cos (#-^ + (,> ) 

- ( 2 S5-6sin(/5 2 ) - ( 2 w5_ 6 cos(/3 2 ) (8) 

where a (=dw/dt) is the angular acceleration of the flywheel; S 2-3 
(= dw 2 - 3 /dt) is the angular acceleration of link 2—3; and es_6 
(=d(i >5 6/dt) is the angular acceleration of link 5—6. 


2.3. Piston (1) and its connecting rod (link 2—3) 

The dynamic simulation of all the parts of the mechanism is 
analyzed in this section. According to Fig. 2, for link 2—3, one has 

/isin(^) = RpSin(| —0) (9) 

Differentiating Eq. (9) with respect to time yields 

IlCOS(fr)^ = -J?pWCOs(|-0) (10) 

The angular velocity of the connecting rod 2—3 is then obtained 
as 

W 2 -3 = —Rpwcos(|— 0)//icos(0 1 ) (11) 

Angular acceleration of link 2—3 is obtained by further differ¬ 
entiating Eq. (11) with respect to time as: 

£ 2 _ 3 = ~^cos(/y[acos(|-0) +w 2 sin (|-0) 

-ww 2 - 3 tan(0i)cos(!-0)] (12) 

As the above kinematic equations of the piston and rod 2—3 are 
derived, the dynamic equations are carried out with respect to the 
moving frame (n—t plane), and the torque balance equation can be 
written in terms of the crankpin axis as shown in Fig. 2. That is, 

Force balance in n-direction: 


f 2n - WjeosO?!) +F 3n -F lx sin(^)-fiyCOS(^) = 0 (13) 

Force balance in t-direction: 


F 2t - W^in^)—F 3t + F ly sin(^)-F lx cos(/3i) = 0 (14) 

Torque balance with respect to joint 3: 


-7’p+^ 2 p/i+jW 1 / 1 sin(ft)+F lx / 1 sin(ft)-Vi C0S (^) = ° ( 15 ) 


Fp — 7n*2 (P c — F a tm) 

(16) 

W P = rn p g 

(17) 

F\y = — Fp — Wp + ITlpClp 

(18) 

F 2 n = 2 1 m 2-3^2-3 

(19) 

c 1, 

Fit = 2*l m 2-3 w 2-3 

(20) 

Wi = m 2 -3 (g - Qp) 

(21) 



where 


Fig. 3. Free-body diagram of displacer and its connecting rod. 
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Tp = l2-3 e 2-3 (22) 

Note that F p is the periodic gas force caused by the pressure in 
the compression chamber which can be determined from the 
thermodynamic model [23], W p is the weight of the piston 
assembly, and VIA, is the weight of the rod 2—3 plus the inertia 
effects caused by the piston’s reciprocating movement. Equations 
(13)—(15) could be solved for the three unknown forces: Fi x , F 3n , 
and F 3t . 

The tangential force exerted by the piston on the flywheel can 
then be written as 

Ftin = -F 3 nSin(|-0 + 0i) +F 3t cos(|-0 + 0i) (23) 

2.2. Displacer (4) and its connecting rod (link 5—6) 

The free-body diagram of the displacer and its connecting rod 
5—6 are plotted in Fig. 3. For the displacer and its connecting rod, 
one has 


Z 2 cos(/ 3 2 ) C ^ = R d wc os(0-|+ (25) 

The angular velocity of the connecting rod can be given as 
W 5-6 = R d wcos(0-~ + $)// 2 cos(|3 2 ) (26) 

The angular acceleration of the connecting rod is 

£ 5—6 = os(/J 2 ) [acos (d - ^ - w 2 sin (d - ^ 

-ww 5 _ 6 tan(0 2 )cos(0-|+$)] (27) 

Then, the force and torque balance equations for the displacer 
and its connecting rod are written with respect to crankpin axis as 
shown in Fig. 3 as follows: 

Force balance in n-direction: 


/ 2 sin(/? 2 ) = R d sin(0-^+$) 

Differentiating Eq. (24) with respect to time yields 


(24) 


-N 2n - W 2 cos(fe) + JV 3n + N u sin(fo) - N, y cos(/? 2 ) = 0 (28) 


Force balance in t-direction: 



-JV 2t + W 2 sin(fe) - N 3t + N ly sin(fo) + N lx cos(/3 2 ) = 0 (29) 


Torque balance with respect to joint 3: 


-I'd + ^N 2f ( 2 - ilV 2 / 2 sin(fo) - JV ly / 2 sin(ft,) - N lx ( 2 cos(/3 2 ) 


= 0 

(30) 

where 


F d = 7tr?(P e - Pc) 

(31) 

W d = m d g 

(32) 

Nly = -Fd - W d + m d a d 

(33) 


Table 1 

Geometrical variables of the baseline case. 


Variable 

Value 

Variable 

Value 

R p (m) 

0.0095 

a>i (rpm) 

1500 

R d (m) 

0.0067 

Pi (kPa) 

101.3 

h (m) 

0.079 

Th (K) 

1300 

h (m) 

0.055 

T C (K) 

300 

fp(m) 

0.0389 

<P (degree) 

70 

Id (m) 

0.045 

dr (s) 

0.000001 

r, (m) 

0.0095 

f'eff 

0.75 

r 2 (m) 

0.0099 

Re,i (°C/kW) 

40118.312 

C(m) 

0.0004 

R cj ("C/kW) 

40689 

let (m) 

0.055 

Gas medium 

Air 

ft(m) 

0.1724 



m 2 -3 U<g) 

0.002187 



m 5 -6 (kg) 

0.001563 



m p (kg) 

0.014013 



m d (kg) 

0.013292 



/(kg-m/s 2 ) 

0.0008 




Fig. 4. Flow chart of the computational process. 
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Table 2 

Operating conditions and the results of steady-state rotational speed, output power, and thermal efficiency of studied cases. 


Case 

u>t (rpm) 

P (kPa) 

Th (K) 

<P (degree) 

£ eff 

/(kg-m 2 ) 

C(m) 

R p (m) 

id (m) 

(L> ave (rpm) 

W(W) 

V 

1 

1500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

899.63 

7.9089 

0.3580 

2-1 

1000 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

899.10 

7.9087 

0.3581 

2-2 

500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

898.18 

7.9093 

0.3581 

3-1 

1500 

101.3 

1000 

70 

0.75 

0.0001985 

0.0004 

0.0067 

0.045 

887.01 

7.7958 

0.3567 

3-2 

1500 

101.3 

1000 

70 

0.75 

0.00008 

0.0004 

0.0067 

0.045 

849.27 

7.4820 

0.3521 

3-3 

1500 

101.3 

1000 

70 

0.75 

0.004 

0.0004 

0.0067 

0.045 

1210.24 

7.7870 

0.3324 

3-4 

1500 

101.3 

1000 

70 

0.75 

0.008 

0.0004 

0.0067 

0.045 

1335.44 

7.6209 

0.3203 

4-1 

1500 

202.6 

1000 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

996.34 

8.9492 

0.3548 

4-2 

1500 

303.9 

1000 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

1026.66 

9.2891 

0.3532 

5-1 

1500 

101.3 

800 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

631.76 

5.2236 

0.3624 

5-2 

1500 

101.3 

1300 

70 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

1247.30 

11.8237 

0.3512 

6-1 

1500 

101.3 

1000 

40 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

933.56 

6.6530 

0.2938 

6-2 

1500 

101.3 

1000 

50 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

966.56 

7.5149 

0.3289 

6-3 

1500 

101.3 

1000 

60 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

944.66 

7.8983 

0.3496 

6-4 

1500 

101.3 

1000 

80 

0.75 

0.0008 

0.0004 

0.0067 

0.045 

845.88 

7.6635 

0.3569 

7-1 

1500 

101.3 

1000 

70 

0.76 

0.0008 

0.0003 

0.0067 

0.045 

843.41 

7.2848 

0.3404 

7-2 

1500 

101.3 

1000 

70 

0.72 

0.0008 

0.0005 

0.0067 

0.045 

818.28 

7.0795 

0.3174 

7-3 

1500 

101.3 

1000 

70 

0.69 

0.0008 

0.0006 

0.0067 

0.045 

726.15 

6.1572 

0.2761 

8-1 

1500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0070 

0.043 

906.85 

7.9903 

0.3565 

8-2 

1500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0070 

0.047 

960.88 

8.5622 

0.3915 

9-1 

1500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0070 

0.045 

935.27 

8.2887 

0.3743 

9-2 

1500 

101.3 

1000 

70 

0.75 

0.0008 

0.0004 

0.0064 

0.045 

862.90 

7.5258 

0.3419 


N 2n = 2*2 m 5-6 w 5-6 

(34) 

N 2t = 2*2 m 5-6 £ 5-6 

(35) 

W 2 = m 5 _ 6 (g~a d ) 

(36) 

Td = h-6 e 5-6 

(37) 


In the above equations, Fd is the gas force caused by the pres¬ 
sure difference between expansion and compression chambers, 
which is calculated from the thermodynamic model [23], Wd is the 
weight of the displacer assembly, Td is the torque on the center of 
gravity of the rod, and W 2 is the weight of the rod 5—6 plus the 
inertia force caused by the displacer’s reciprocating movement. 
Equations (28)—(30) are solved for three unknown forces: Ni x , N 3 „, 
Ni t . The tangential force exerted by the displacer on the flywheel 
can then be calculated as 

Ntin = W 3n sin(6»-|+0+fe)-W3tCOs(6»-|+<P+fo) (38) 

2.3. Equation of rotational motion of flywheel 

It is assumed that the flywheel dominates the moment of inertia 
of the engine; therefore, the rotational speed of the engine can be 
determined by dealing with the rotational motion of the flywheel. 
The tangential forces, F t i n and N t j n , calculated by Eqs. (23) and (38), 
respectively, are introduced to determine the gas pressure torques 
exerted on the flywheel as: 

Tp = Ftinflp (39) 

T d = N tin R d (40) 

On the other hand, the friction torque in the engine is also taken 
into account in this study. Note that a number of existing studies 


[26,27] dealt with the friction model of the internal combustion 
engines. However, these friction models cannot be used for the 
beta-type Stirling engines whose structure is different from the 
internal combustion engines. In a recent report, Cheng and Yang 
[28] proposed an empirical correlation of the friction work loss in 
terms of the angular speed of a beta-type Stirling engine. According 
to the observation of [28], the friction torque is assumed to be 
proportional to the mechanical loss of the engine as 

Tf r = c fr w (41) 

where as a test case, the friction factor (cf r ) can be derived from the 
empirical correlation for the engine designed by Cheng and Yang 
[28] 

c fr = 2.98 X 10 8 + 1.02 X lO- 4 ^- 1 (42) 

It is noticed that the coefficients of Eq. (42) are dependent on 
individual Stirling engine in use. A universal correlation for the 
friction factor suitable for all different types of Stirling engines is 
not realistic. For other type of Stirling engines, the form of friction 
model could even be altered dramatically. 

The conservation of angular momentum is applied to yield the 
equation of rotational motion as 

Idco/dt = T p + Tj — Tf r — T] oac |. (43) 

where / is the total mass moment of inertia of the engine’s 
flywheel and Ti oa( j represents the external load of the engine. Eq. 
(43) can be integrated numerically with respect to time by means 
of the fourth-order Runge—Kutta method [29] to yield the tran¬ 
sient variation in the rotational speed w. Note that in the present 
dynamic model, the effects of friction and components' weights 
and mass inertia are included. Also, it is important to mention 
that the external load Toad is treated as a variable in the present 
model which can be specified in accordance with the practical 
loading condition. Here, the value of Toad is assigned to be zero 
so as to investigate the balance between the gas pressures and 
the friction force. 
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t(sec) 

Fig. 5. Transient variations in pressures and temperatures in the expansion and the 
compression chambers, for baseline case (case 1). 

3. Results and discussion 

In the numerical simulation process, the engine is started from 
an initial rotational speed. The torques exerted by the flywheel of 
the engine at any time instant caused by the gas pressures in the 
chambers, the mass inertia, the friction force, and the external load 
are calculated. The instantaneous rotation speed of the engine is 
then determined by integration of the equation of rotational 
motion, Eq. (43), with respect to time. The obtained instantaneous 
rotation speed of the engine is introduced into the thermodynamic 
model to determine the instantaneous variations in pressure and 
other thermodynamic properties of the gas inside the chambers. 
The gas pressures are used to calculate the gas pressure torques, 
and the computation marches to the next time step. The transient 
variations in gas properties inside the engine chambers and the 
dynamic behavior of the engine mechanism are handled simulta¬ 
neously via the coupling of the thermodynamic and dynamic 
models. Fig. 4 shows the flow chart of computation process. The 



Fig. 6. Instantaneous and cycle-average rotational speeds from different initial speeds 
(cases 1, 2—1, and 2-2). 

geometry specification and operating conditions of a baseline case 
are listed in Table 1. In the baseline case, initially the gas medium, 
air, is assumed to be 300 K in both chambers and the charged 
pressure is 101.3 kPa. The temperature of the wall of expansion 
chamber is maintained at T H from the beginning of the simulation 
since the hot-start simulation is of major concern here. The value 
of T h is changed in the parametric study whereas T L is fixed at 
and 300 K. In addition to the baseline case, there are 20 more 
cases considered in a comprehensive parametric study. The oper¬ 
ating conditions of these studied cases are given and numbered in 
Table 2, and the obtained results of the rotational speed, output 
work, and thermal efficiency are provided. 

Fig. 5 shows the transient variations in pressures and temper¬ 
atures of air in the expansion and compression chambers (P e . Pc, T e , 
and T c ). As shown in Fig. 5(a), the pressures in the expansion and 
the compression chambers are very close. It is reasonable because 
the engine is small so that the pressure within the cylinder of the 
engine is almost equal. When the regenerative channel gap is 
reduced, a larger pressure drop between the chambers may be 



t (sec) 

Fig. 7. Effects of moment inertia of flywheel on instantaneous and cycle-average 
rotational speeds (cases 1, 3—1, 3-2, 3—3, and 3—4). 








































































































































































C.-H. Cheng, Y.-J. Yu / Renewable Energy 36 (2011) 714-725 


721 



Fig. 8. Transient variations in cycle-average output power with different moment 
inertia of flywheel (cases 1, 3-1, 3-2, 3-3, and 3-4). 


expected. In addition, the oscillating variation of the pressures 
reflects the periodic movement of the piston and displacer. The 
average pressure reaches a value of 200 kPa in about 0.5 s from 
start. The pressures inside the chambers lead to the gas force acting 
on the piston, displacer, and hence flywheel. In the second plot of 
this figure, the transient temperature variation in both chambers is 
depicted. The temperatures in the two chambers (T e and T c ) are 
oscillatory and increased from 300 I< to approximately 700 K and 
600 K in the expansion and the compression chambers, respec¬ 
tively. Note that the difference between T H and T e in the steady 
regime is around 300 K and the difference between T c and Tl is 
nearly of the same magnitude. The temperature differences (Th - T e 
and T c - Tl) are strongly related to the heats absorbed and rejected 
per cycle. 

Fig. 6 shows the transient variation in the rotational speed of the 
engine starting from different initial rotational speeds ((a, = 500, 
1000, 1500 rpm). The cycle-average rotational speed (&j ave ) is also 
calculated and is illustrated with black curves. It is observed that 
the initial rotational speed indeed has subtle influence on rota¬ 
tional speed in the start-up period; however, after the start-up 
period all different initial rotational speeds lead to the same final 
cycle-average rotational speed in the steady-state regime. This is 
reasonable because from Eq. (43) as in the steady-state regime, the 
time derivatives of the cycle-average quantities vanish, that is 
d(t> ave /df= 0. As a result, the steady-state cycle-average rotational 
speed finally reached should be 

Wave = (7p, av e + ^d.ave — ^load,ave)/Cfr (44) 

It implies that in the steady-state regime the cycle-average 
rotational speed (w s ) is independent of the initial rotational speed. 
However, it is observed that only an initial rotational speed within 
a certain range can start the engine and make the engine reach the 
steady-state regime finally. The engine is brought to a stop if an 
improper initial rotational speed outside the range is chosen. 

Fig. 7 conveys the effects of moment inertia of flywheel on the 
instantaneous and the cycle-average rotational speeds, for cases 1, 
3—1, 3—2, 3—3, and 3—4. In this figure, the magnitude of the 
moment inertia of flywheel is varied from 0.00008 to 0.008 kg-m2. 
It is found that the oscillating amplitude of the instantaneous 




Fig. 9. Effects of initial charged pressure on cycle-average rotational speed and output 
power (cases 1, 4-1, and 4-2). 

rotational speed is decreased by an increase of the moment of 
inertia of the flywheel. Meanwhile, for the case with a flywheel of 
larger moment inertia, it takes longer time for the engine to reach 
the steady-state regime. However, as the steady-state regime is 
reached, the cycle-average rotational speed ((<j s ) approaches the 
same value regardless of the magnitude of moment inertia of 
flywheel. 

In Fig. 8, the dependence of engine’s cycle-average output 
power (W) on the moment of inertia of the flywheel is depicted for 
the same cases. In this figure, it is observed that the cycle-average 
output power in the steady-state regime attains a maximum of 
7.91 W (and the thermal efficiency reaches a maximum of 0.358 as 
given in Table 2) when the value of I is equal to 0.0008 kg-m 2 . This 
obviously implies that for any particular engine design, there exists 
an optimal value for the moment of inertia of the flywheel that 
should be determined in order to maximize the output power and 
the thermal efficiency. 

Fig. 9 displays the effects of initial charged pressure on the 
cycle-average rotational speed and the cycle-average output 
power of the engine, for cases 1, 4—1, and 4—2. It is found that 
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both the cycle-average rotational speed and the cycle-average 
output power increase with the initial charged pressure. A 
higher initial charged pressure means a larger mass of air being 
charged in the cylinder of the engine. With more gas medium 
in the engine, the output power is expected to increase. The 
initial charged pressure is assigned to be 101.3, 202.6, and 
303.9 kPa here. In Fig. 9 and Table 2, it is found that the power 
out achieves a maximum of 9.2891 W at P = 303.9 kPa. However, 
based on the thermal efficiency data shown in Table 2 for cases 
1, 4—1 and 4—2, it is observed that an increase in initial charged 
pressure may lead to a decrease in the thermal efficiency. The 
reason for the decrease in the thermal efficiency accompanying 
the increase in the charged pressure is attributed to the higher 
heat input required for the engines containing more mass of air. 
An increase in mass of air elevates the thermal inertia of the 
gas medium; therefore, higher heat input is required to heat up 
the gas medium. In consequence, elevating the charged pressure 
of the engine could increase the engine’s output power but 
might reduce the thermal efficiency of the engine. 




Fig. 10. Effects of heat source temperature on cycle-average rotational speed and 
output power (cases 1, 5—1, and 5—2). 


Fig. 10 shows the heat source temperature effects on the 
performance of the engine. The heat source temperature is set at 
800 K, 1000 K, or 1300 K. It is observed that the rotational speed 
and the output power are elevated by increasing the heat source 
temperature. However, the value of Th is strictly restricted by the 
thermal expansion of the materials of the heating head. Since the 
temperatures of thermal reservoirs are given, the heat transfer to 
the chamber can be calculated in terms of the thermal resistance 
and the temperature difference between the fluid temperature and 
the heat source temperature. A higher heat source temperature 
means a higher heat transfer rate input to the expansion chamber. 

Fig. 11 shows the effects of the phase angle between the 
displacements of the displacer and the piston on the cycle-average 
rotational speed and output power. The phase angle of the baseline 
case (case 1) is assigned to be 70 degrees. The simulation is made for 
four more phase angles, 50, 60, 80, and 82 degrees (cases 6—1, 6—2, 
6—3, and 6—4) for comparison. This figure indicates that the optimal 
values of the phase angle leading to highest rotational speed or 
output power exist. According to the data provided in Fig. 11 and 




Fig. 11. Effects of phase angle on cycle-average rotational speed and output power 
(cases 1, 6—1, 6-2, 6-3, and 6-4). 
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Table 2, it is found that the case at 0 = 50 degrees leads to the 
highest cycle-average rotational speed in the steady-state regime, 
whereas the case at 0 = 70 degrees leads to the highest cycle- 
average output power and thermal efficiency. Note that the effects 
of the phase angle are rather involved. Typically, for the small-scale 
engines with efficient heat transfer between the gas medium and 
the walls of the expansion and the compression chambers, the 
optimal phase angle is set to be less than 90 degrees since the time 
required for heat transfer is shorter. On the contrary, for the larger- 
scale engine with relatively poor heat transfer, the optimal phase 
angle is in general higher than 90 degrees since the gas medium 
takes longer time for sufficient heat transfer. An adjustment in the 
phase angle should be made for different purposes. 

Next, the effects of engine’s geometrical parameters are inves¬ 
tigated. The selected parameters are: regenerative gap size, length 
of the displacer, and the piston stroke. 

Fig. 12 illustrates the subtle influence of the regenerative gap 
size. The regenerative gap size is varied from 0.0003 to 0.0006 m, 
corresponding to cases 1, 7—1, 7—2, and 7—3. As the gap size is 
decreased, significant increase in the shear friction loss is expected. 




Fig. 12. Effects of regenerative gap size on cycle-average rotational speed and output 
power (cases 1, 7—1, 7-2, and 7-3). 


As the gap size is increased, the regenerative effectiveness of the 
regenerative channels will be reduced. Therefore, there should be 
an optimal value for the gap size that leads to a peak value of the 
output power. For the particular case, the gap size of 0.0004 m 
exhibits the highest cycle-average rotational speed and the largest 
output power in the steady-state regime. 

Fig. 13 shows the effects of the length of the displacer on the 
engine’s rotational speed and output power. The rotational speed 
and the output power appear to increase with the displacer’s length 
monotonically. As the stroke of the displacer is fixed, an increase in 
the length of the displacer means a reduction in the dead-zone 
space and an increase in the compression ratio. Therefore, the 
performance of the engine could be improved by increasing the 
displacer’s length as long as the displacer would not collide with 
the walls of the cylinder and the piston. 

Fig. 14 conveys the effects of piston stroke on the cycle-average 
rotational speed and output power, for cases 1, 9—1, and 9—2. The 
stroke of the piston is actually two times R p . Therefore, in this 
figure R p is used as the parameter varied from 0.0064 to 0.0070 m. 
It is found that both the cycle-average rotational speed and the 



Fig. 13. Effects of displacer length on cycle-average rotational speed and output power 
(cases 1, 8-1, 8—2). 
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cycle-average output power increase with R p . This is probably 
because a larger piston stroke leads to a higher compression ratio, 
which typically results in higher output power and thermal 
efficiency. 

4. Concluding remarks 

In this study, dynamic simulation of a beta-type Stirling engine 
with cam-drive mechanism has been attempted. A dynamic model 
associated with the cam-drive mechanism has been developed and 
combined with the thermodynamic model [23] to predict the 
transient behavior of the engine in the hot-start period. Based on 
the numerical simulation results, the following conclusions are 
drawn: 

(1) It is observed that the initial rotational speed and the moment 
inertia of flywheel exhibit subtle influence on rotational speed 
in the start-up period; however, after the period all different 



(b) 

Fig. 14. Effects of piston stroke on cycle-average rotational speed and output power 
(cases 1, 9—1, and 9—2). 


initial rotational speeds and moment inertia of flywheel lead to 
the same final cycle-average rotational speed in the steady- 
state regime. In addition, for the particular cases, it is observed 
that the cycle-average output power and the thermal efficiency 
attain their respective maximums when the value of moment 
inertia of flywheel is equal to 0.0008 kg-m 2 . 

(2) Elevating the charged pressure of the engine could increase the 
engine’s output power but might reduce the thermal efficiency 
of the engine. 

(3) Heat input to the expansion chamber increases with the heat 
source temperature. Therefore, it is observed that the rotational 
speed and the output power are elevated by increasing the heat 
source temperature. However, the value of Th should be strictly 
restricted by the thermal expansion of the materials of the 
heating head. 

(4) Typically, for the small-scale engines with efficient heat 
transfer between the gas medium and the walls of the expan¬ 
sion and the compression chambers, the optimal phase angle is 
set to be less than 90 degrees since the time required for heat 
transfer is shorter. In the simulation, totally five angles, 50, 60, 
70, 80, and 82 degrees are taken into consideration in 
comparison. It is found that the case at <P = 50 degrees leads to 
the highest cycle-average rotational speed in the steady-state 
regime, whereas the case at $ = 70 degrees leads to the highest 
cycle-average output power and thermal efficiency. 

(5) The effects of engine’s geometrical parameters, including 
regenerative gap size, length of the displacer, and the piston 
stroke, are also investigated for the baseline case. Results show 
that the gap size of 0.0004 m exhibits the highest cycle-average 
rotational speed and the largest output power in the steady- 
state regime. In addition, since an increase in the length of the 
displacer or in the piston stroke leads to a reduction in the dead- 
zone space and an increase in the compression ratio, the 
performance of the engine could be improved by increasing the 
displacer’s length or the piston stroke. 
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